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ABSTRACT 



We investigate the viability of the magnctorotational instabihty (MRI) in X- 
ray ionized viscous accretion disks around both solar-typc stars and very low 
mass stars. In particular, we determine the disk regions where the MRI can be 
shut off either by Ohmic resistivity (the so-called Dead and Undead Zones) or by 
ampipolar diffusion (a region we term the Zombie Zone). We consider 2 stellar 
masses: = 0.7 Mq and O.IMq. In each case, we assume that: the disk surface 
density profile is that of a scaled Minimum Mass Solar Nebula, with Mdisk/M^ — 
0.01 as suggested by current data; disk ionisation is driven primarily by stellar 
X-rays, complemented by cosmic rays and radionuclides; and the stellar X-ray 
luminosity scales with bolometric luminosity as Lx/L^, ^ 10~^'^, as observed. 
Ionization rates are calculated with the MOCCASIN Monte Carlo X-ray transport 
code, and ionisation balance determined using a simplified chemical network, 
including well- mixed 0.1 /im grains at various levels of depletion. We find that 
(1) ambipolar diffusion is the primary factor controlling MRI activity in disks 
around both solar-type and very low mass classical T Tauri stars. Assuming that 
the MRI yields the maximum possible field strength at each radius, we further find 
that: (2) the MRI-active layer constitutes only ~5-10% of the total disk mass; 
(3) the accretion rate (M) varies radially in both magnitude and sign (inward or 
outward), implying time- variable accretion as well as the creation of disk gaps 
and overdensities, with consequences for planet formation and migration; (4) 
achieving the empirical accretion rates in solar-type and very low mass stars 
requires a depletion of well-mixed small grains (via grain growth and / or settling) 
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by a factor of 10-1000 relative to the standard dust-to-gas mass ratio of 10~^; 
and (5) the current non-detection of polarized emission from field-aligned grains 
in the outer disk regions is consistent with active MRI at those radii. 



Introduction 



The inward migration of material in protoplanetary disks requires a source of anoma- 
lous viscosity; microscopic gas viscosity is far too inefficient to explain the observed accretion 
rates and disk lifetimes. The best current candidate is turbulent viscosity generated by the 
magneto-rotational instability (MRI), which arises due to Keplerian shearing of weak mag- 



netic fields embedded in a partially ionized disk (Balbus & Hawley 1991). Assuming ioniza- 
tion due to cosmic rays, Gammie (1996) argued that MRI-driven accretion in protoplanetary 
disks is likely to be "layered" : near the mid-plane in the inner disk, where column densities 
are too high for cosmic rays to penetrate, Ohmic resistivity quenches the MRI and forms 
a Dead Zone, while active MRI proceeds in the surface layers of the inner disk as well as 
throughout the outer disk, where low densities permit sufficiently high fractional ionization 



for good gas-field coupling. Igea & Glassgold (1999, hereafter IG99) subsequently showed 



that ionization by stellar X-rays overwhelms that by cosmic rays, but this does not change 
the basic paradigm of an Ohmic-dominated Dead Zone hemmed in by active disk regions. 

It has become increasingly clear, though, that this picture is overly simplified, because 
Hall and ambipolar diffusivities cannot be ignored under the conditions prevalent in these 
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). Specifically, active MRI in this 



possible where the density is low enough to allow adequate fractional ionization (i.e., outside 
the Dead Zone); however, it is precisely at these low densities that either electrons, or both 
electrons and ions, can become decoupled from the neutrals, promoting Hall diffusion in 
the first instance and ambipolar in the second. The non-dissipative nature of the Hall term 
means that it may enhance, suppress or leave unaffected the MRI efficiency, depending on 



the precise environment and magnetic field geometry (Wardle 1999; Sano & Stone 2002 
Wardle fc Salmeron||2012 ). The ambipolar term, though, is dissipative, and can smother the 
MRI by decoupling the field entirely from the predominantly neutral gas. Dust grains do 
not change these conclusions, but only add to the dangers facing the MRI, by soaking up 
electrons and facilitating ion-electron recombinations on their surfaces, thus greatly reducing 
the ionization fraction and increasing the resistivities (e.g., Sano et aL||2000 Ilgner & Nelson 



2006 hereafter IN06). 



In a seminal analysis. 



Bai & Stone 



(2011) have investigated in detail the effects of 



3 



ambipolar diffusion on the non-linear evolution of the MRI, through local 3D shearing-box 
simulations. Their calculations are in the "strong coupling" limit (wherein the ion density 
is negligible compared to the neutral density, and the electron recombination time is much 
smaller than the orbital period), and explicitly account for the non-conservation of ion densi- 
ties (which depend instead on the chemical equilibrium). Their analysis yields a quantitative 
general criterion for evaluating whether or not ambipolar diffusion can extinguish the MRI. 



Following this, Bai (2011, hereafter Bll) showed that the "strong coupling" limit does 
invariably hold in protoplanetary disks. Applying the above criterion to X-ray irradiated 
disks around solar-type stars, he further finds that ambipolar diffusion strongly determines 
the MRI efficiency in such disks, and, in the presence of small (0.1 //m) grains, drives the 
attendant accretion rate below the observed rates in these stars by an order of magnitude 



or more. In an elegant and complementary semi-analytic study, Perez-Becker fc Chiang 



(2011a, hereafter PCll) find an even graver suppression of the MRI with polycyclic aromatic 
hydrocarbon (PAH) grains. These somber implications for the viability of the MRI are 



alleviated somewhat by Perez-Becker & Chiang (2011b), who contend that stellar FUV 



photons produce much more ionized surface plasmas than X-rays, yielding efficient MRI and 
adequate accretion rates beyond ~10 AU in solar-type disks, in spite of PAHs. Even in this 
case, however, the inferred accretion in the inner disk remains much weaker than observed. 

Our goals in this paper are twofold. First, we wish to extend the above studies of X-ray 
irradiated disks to very low mass stars {^0.1 Mq, i.e., mid-M spectral types). It is only 
by expanding the parameter space in this fashion that we can hope to fully understand the 
feasibility and importance of the MRI for driving viscous accretion. For instance, it has been 
suggested that the disks girdling such stars may be far more active than their more massive 
counterparts around solar-type stars, because their lower column densities enable ionizing 



photons to penetrate much deeper (e.g., Hartmann et al. 2006). In that case, perhaps the 



problems, alluded to above, that plague the MRI in more massive disks are not as severe for 
very low mass stars? However, as discussed, lower densities also enhance ambipolar diffusion, 
so it is a priori not at all clear if this hypothesis is valid (as we will see, it is not). Moreover, 
there are now burgeoning observations of planets down to terrestrial masses orbiting low- 
mass M dwarfs (e.g., Bonfils et al.||2011 ), and understanding their formation and migration is 
inextricably linked to deciphering the viscosity and accretion properties of their natal disks. 

Second, we wish to implement certain important improvements over the previous anal- 
yses. In particular, we undertake a more sophisticated treatment of the X-ray ionization, 
including Compton scattering, with the 3-D Monte Carlo photoionization and dust radiative 
transfer code moccasin. We also derive accretion rates via a more generalized equation, 
and show in the process that radial gradients have some important implications (e.g., time- 
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dependent accretion rates, and the formation of disk gaps and overdensities). Furthermore, 
we examine the effects of various levels of grain depletion, intermediate between the limiting 
cases of a standard ISM dust abundance and no dust at all. Finally, for both solar-type 
and very low mass stars, we compute the detailed chemical, ionization and MRI properties 
over the entire disk from 0.1 to 100 AU, instead of at a selected few radial locations. Save 
these changes (and some differences in the chemical network, discussed at the appropriate 
junctures), our study is analogous to that of Bll. 

The layout of the paper is as follows. We outline the criteria for active MRI in §2, and 
discuss the conditions for determining the relevant field strengths and accretion rates in §3. 
In §4, we describe the adopted disk physical properties - the disk structure (§4.1), chemical 
and ionization balance (§4.2), and resulting resistivities (§4.3) - used to evaluate the MRI 
criteria given in §2. The specific stellar cases we examine are summarized in §5. Our results 
are presented in §6, and conclusions and directions for future work laid out in §7. 

2. Conditions for Magnetically Driven Disk Accretion 
2.1. Dead and Undead Zones 



Magnetic stresses can transport angular momentum in the disk, and thus drive disk 
accretion, only if the gas is sufficiently ionized to couple to the field, i.e., if gas motions can 
generate magnetic disturbances faster than these can diffuse away due to the finite resistivity 
f]. In a Keplerian disk, it is fundamentally the orbital shear that creates the stresses, so the 
relevant minimum timescale for field generation is the orbital period, ~ For large-scale 
toroidal fields produced by such shearing of a weak seed radial field, the characteristic length- 
scale for field variations is of order the radial distance to the star, ~ vk/^, where vk is the 
Keplerian velocity. The dissipation timescale is thus ~ {vk/^)'^ /i]o, if Ohmic resistivity r/o 
provides the primary diffusion mechanism. The condition for field regeneration to prevail 



over dissipation is then given by the Reynolds number criterion (Turner & Sano 2008): 



with the lower limit chosen to give toroidal fields 10 times stronger than the initial seed. 

For local tangled fields generated by MRI-driven turbulence, the height of the thin 
disk sets an upper limit on the wavelength of vertical MRI modes, and thereby on the 
dissipation timescale, so it is the vertical direction we must consider. For a vertical mode with 
wavenumber k, the Ohmic dissipation rate is ~ k'^rjo while the growth rate is kvj^, where vj^ 
is the local Alfven velocity. The maximum growth rate is again the orbital frequency ~ Q; 
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the fastest growing mode thus has a wavenumber k = Q/vj^. Demanding that its growth 
rate exceed its dissipation rate then yields the Elsasser number criterion: 



A 



A 



> 1 



(2) 



whether the net background field is vertical, toroidal or zero (Sano & Stone 2002). As 



noted by Turner et al. (2007), the above Elsasser number criterion captures MRI activity 
more accurately than an analogous version with replaced by the sound speed Cs (e.g.. 



as adopted by Fleming et al. 2000 Wardle 2007). The specification of the Alfven velocity 
v_A[r] = B[r]/ ^yA^^p\r] is described in §3. For now, note two consequences of the requirement 
(see §2.2 below) that the gas pressure dominate over the magnetic pressure in the disk. First, 
this implies Vjs,/cs <^ 1, which guarantees that the wavelength ~ Vjs,/Vt of the fastest growing 
mode near the midplane is indeed much less than the disk scale height zh ~ Cs/f2. 

Second, combining the requirement vj^jcg ^ 1 with the thin disk condition vk ^ Cs 
implies vk ^ vj^. Hence, of the two criteria (1) and (2), the first is the less restrictive, in 
that it allows field generation at much higher resistivities than the second. Thus the two 
divide the disk into 3 parts. In the active zone, both (1) and (2) are satisfied, and efficient 
magnetically-driven angular momentum transport can occur in spite of Ohmic resistivity. 
In the Undead Zone, the MRI is quenched, but orbital shear can still generate large-scale 
toroidal fields out of radial ones siphoned off from neighbouring active layers (if any exist). In 
the Dead Zone, both processes are quenched by Ohmic dissipation, the gas is almost entirely 
decoupled from the field, and magnetic torques cannot drive angular momentum transport. 



2.2. Ambipolar Diffusion and the Zombie Zone 

The criteria in equations (1) and (2) have been formulated in terms of the Ohmic 
resistivity. However, Hall or ambipolar diffusivities can be expected to dominate over Ohmic 
in various parts of the disk. What are the criteria for MRI under such conditions? 

Equation (2) simply expresses the condition that, for robust MRI in the Ohmic domain, 
the ratio of the inductive to the Ohmic resistive term in the induction equation must exceed 
unity. When Hall or ambipolar diffusivity is the major resistive term instead, the analogous 
condition is naively obtained by replacing r]o in equation (2) by \'r]H\ (Hall) or r]A (ambipolar) 



(see §4.3). However, local 3-D MHD simulations by Sano & Stone (2002) show that, while 
the Hall term can affect the saturation level of the Maxwell stress by a factor of a few, 
it hardly changes the critical condition for strong MRI: this remains the Ohmic Elsasser 
criterion equation (2), regardless of the size of the Hall term or the strength and geometry of 
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the field. Hence we use equation (2) in the Hall domain as well (i.e., calculate the Elsasser 



number using 770 even when \riH\ dominates, thereby ignoring Hall effects altogether). Wardle 



& Salmeron (2012) show that the aforementioned simulations do not probe all the pertinent 



parameter space, and perform a linear analysis that suggests Hall diffusion may have a 
powerful effect on the MRI in a previously unexplored parameter regime that is nevertheless 
very germane to protoplanetary disks. We discuss their results briefly in §7; however, in the 
absence of non- linear simulations here, we neglect this added complication for the present 
analysis. 



In the ambipolar regime, on the other hand, Wardle (1999) argues that it is appropriate 



to substitute r]A for r/o in equation (2), leading to the ambipolar Elsasser number condition 



Am = 



> 1 



(3) 



for the MRI to exist here. Moreover, if electrons and ions are the only charged species, this 
reduces to 7jPi/fi > 1 (see §4.3), where 7^ is the neutral-ion collisional drag coefficient and pi 
the ion density; this implies that the MRI can occur under these circumstances if a neutral 
particle collides with at least one ion on average per orbital period. This condition has been 



invoked by several authors (e.g.. Turner et al. 2010) in ambipolar-dominated disk regions. 



Until recently, the only existing simulations of MRI with ambipolar diffusion, by Haw 



ley & Stone (1998), indicated that the above criterion was actually too liberal. Using 3-D 



MHD local shearing box simulations, based on an idealised 2-component fluid (ions + neu- 
trals), these authors found that efficient angular momentum transport by the MRI requires 
neutral-ion collisions to be at least a hundred times more frequent, i.e., jipi/Q > 100. For 
'jipi/fl ~ 1, the linear instability is still present, but the nonlinear evolution of the MRI is 
controlled primarily by the ion density, and the total angular momentum transport becomes 
proportional to the ionisation fraction. In other words, when 'jiPi/Q ~ 1, the MRI is alive 
only in the ions, and effectively dead in the neutrals (with the neutrals only providing a drag 



that reduces the angular momentum transport in the ions; Hawley & Stone 1998). 



However, very recent simulations by Bai & Stone (2011), based on a multi-component 
fluid including grains (very similar to the system adopted here; §§4.2, 4.3), suggest that 
the above results are overly pessimistic. With the ratio of the inductive to the ambipolar 
term now given by the full general value v'^/rjA^ (instead of the reduced value •jiPi/fl), 
these authors flnd that the MRI can be sustained at any value of Am, provided the fleld 
is sufficiently weak. Specifically, the MRI can operate as long as the plasma /3-parameter, 
/3 = Pgas/Ps, satisfies 



/3>/3n 



(4a) 
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where the limiting value is given by 



50 



Am 



1.2 



+ 



1/2 



(46) 



with Am given by equation (3) and Pb = B^/Sn. We follow Bll in adopting the above 
condition for MRI in ambipolar-dominated regions. If equation (4) is not satisfied in these 
regions, then the MRI still drives angular momentum transport within the charged species 
that are tied to the field, but is quenched in the neutral bulk of the disk, since the neutrals 
no longer couple to the field through collisions with charged particles. As such, the net 
transport of angular momentum via the MRI becomes negligible. We call this "marginally 
alive, effectively dead" region the Zombie Zone. 

Finally, note that /3min approaches {50 /Am^-"^) for Am < 1, and asymptotes to 1 from 
above as Am — >■ oo. Thus the above condition for sustainable MRI demands that the gas 
pressure dominate over the magnetic pressure in the disk, as stated earlier in §2.1. 

To summarize, our conditions for active MRI are: Ak > 10, A_4 > 1 and /3 > /3min- Note 
that /3min is a function of Am, which itself depends on the field strength through both the 
Alfven velocity and rjA (see §4.3). Similarly, the extent of the Dead Zone also depends on 
the magnetic field via the Alfven velocity. We describe below our method for determining 
the field strength, based on considerations of the accretion rate. 



3. Accretion Rate 

For radial angular momentum transport through magnetic torques, the accretion rate 



(positive in the inward direction) at any radius is given by (Wardle 2007): 

Vk or ^ 

where vk = tVL is the Keplerian velocity, and we have defined 

{—BrB^) = J BrB^dz 



(5a) 



(56) 



with ±s being the upper and lower surfaces of the disk, and 2h the total thickness of the 
MRI-active layers (counting both above and below the disk midplane). Simulations of the 
MRI indicate that B^^^ ~ A^-B^B^) ( Sano et al. 2004 ); we may therefore write: 



M 



1 d 
2rVL dr 



{r'hB^^) 



(6) 
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If M is constant with radius, then the above reduces to (integrating over r, and dropping 
the 'rms' subscript on B for brevity) M = hB^/4Q. This is the formula used by Bll to 
calculate the accretion rate. Notice that this automatically implies that hB^ oc r~^/^ for self- 
consistency with the assumption of a radially invariant M. As we will show, however, there 
is no fundamental reason for hB^ to follow precisely this radial dependence. In particular, 
if this quantity is some arbitrary power-law with radius at a given location, i.e., hB^ oc r~", 
then the general solution to the accretion rate there is M = (2 — n) hB'^/2 For n — 3/2, 
this reduces to the constant M solution above, but for other n the accretion is radially 
variable; moreover, M vanishes or even becomes negative (outward radial flow) for n > 2. 

This brings us to the question of how to set the magnetic field strength. Bll determines 
B by assuming that it is the maximum allowed by the condition in equation (4), /3 = /3min- 
He justifies this on the basis that it allows one to calculate the highest possible accretion 
rate, under the assumption of a radially constant M — hB'^/AO,. As we will show, however, 
M is not constant when B is chosen this way. Nevertheless, wc can still support such a choice 
of B with a more physically motivated justification: that the MRI is maximally efficient, 
generating the strongest field possible that still allows the MRI to operate. 

Our procedure therefore is as follows. Assuming the magnetic field to be vertically 
constant, we loop through a range of field strengths at every location to find the permitted 
-Bmax (corresponding to /3min) at every radius. More precisely, for every value of B (which also 
determines the Alfvcn velocity and the ambipolar resistivity rjj^), wc calculate the thickness 
h of the active layer given our criteria for MRI in §2, and choose the B for which the 
quantity hB^ - the height-integrated stress - is maximized. This is because maximizing 
B alone formally leads to an infinitesimal h, implying M — > 0. However, as we will see, 
HB^ increases in step with B until a maximum value, and then falls off extremely rapidly 
to zero as S is increased further; thus, for all practical purposes, the two maxima coincide, 
i.e., (/i-B^)max occurs essentially at (very shghtly lower than) B ^ Smax- We then insert 
this (/iS^)inax into the general solution, equation [6], to calculate M at every radius. To 
summarize, we set the field strength to maximize the height-integrated stress in the same 
way as Bll does, but our M differ from his since we allow for radial stress gradients. 

To carry out these calculations, we must specify the disk structure, ionization balance 
and chemistry, and the attendant resistivities. These are described below. 
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4. Disk Structure, Ionization Chemistry and Resistivity 
4.1. Disk Structure 



We consider a star of mass M^, radius i?* and bolo metric luminosity L^, surrounded by 
a non-self-gravitating Keplerian viscous accretion disk in vertical hydrostatic equilibrium. 
We assume that the disk is vertically isothermal, with a radial temperature profile set by 
absorbing stellar radiation. We further assume that the disk is spatially thin <^ r, 

where zh[t\ is the vertical pressure scale-height in the disk at a radial distance r from 
the star), equivalent to neglecting all sources of pressure support (e.g., thermal, magnetic) 
compared to rotation in the equation of radial force balance. 

With M^,, L*, r and the mean molecular mass /i of the disk gas normalised as: 



'^*=lMo ' ^* = lZ^' ^=1AU' ^=2:34 



the Keplerian angular rotation velocity becomes 

1/2 

= 2.0 X 10"^ I I s 



r • 



and the radial temperature profile in the disk is 

1/4 



T[r1 = 393|^| K (9) 



r • 



where we have ignored any extra projection factors arising from the relative inclination be- 
tween the disk surface and incident stellar radiation. We note that Hayashi ( 1981[ ) (hereafter 



H81) invokes the same functional form for the temperature, but with a somewhat lower nor- 
malization of 280K based on the assumption that dust is severely depleted in the nebula 
while the gas absorbs only solar UV. Many other recent works examining disk magnetic ac- 



tivity have adopted the latter standard Hayashi formulation (e.g., Bll, Turner et al. 2010). 
Conversely, PCll have employed the temperature structure of a passive flared disk, which 
has nearly the same functional dependence on radius (T oc r~^/^) but a much lower normal- 
isation: ~130K at lAU around a solar mass star. There are a number of uncertainties in 
determining what the disk temperature really is: for example, one most often assumes (as we 
have here) that the disk is passive, with negligible contribution from accretion heating, and 
is also vertically isothermal. Without a rigorous reason for choosing a precise normalisation, 
we have opted here for the most naive estimate. Since the reaction rates (§4.2, Table 1) go as 
the square-root of the temperature, the difference between our normalization and e.g. that 
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of H81 leads to negligible changes in the results, as we explicitly show in our comparisons to 
the results of other groups in §6.1. The sound speed is then 



Cs[f] = 6.1 X 10 



cms 



1.2 X 10' 



cms 



(10) 



so that the disk pressure (and density) scale-height becomes 



Zh\r 



V2 



Cs[r\ 



8.5 X 10 



11 



10 



1/8 



cm 



:iii 



Furthermore, the mass of a disk with a surface density profile S[r] is given by: 

(27rr S[r] dr) 



:i2) 



We adopt a disk outer radius Rout = constant = lOOAU over the range of stellar masses 
investigated here. With insufficient data at present about how this parameter changes with 
stellar mass, this assumption, which sets Rout to be of order that seen in classical T Tauri 
stars, appears reasonable. We also adopt a disk inner radius ~ 5-R*, in agreement with 
magnetospheric truncation models and observations; the precise value of this parameter is 
immaterial to our calculations as long as Rm lies inside the radius where collisions of alkali 
metal atoms with other species take over as the main ionization process, that is, where 
thermal ionization dominates. Finally, we assume that the disk mass varies as oc M^,. 
This relationship, albeit with large scatter (some of which may be age-related) is consistent 
with current data (Scholz et al. 2006 Mohanty et al. 2012). Specifically, we adopt Ma = 
10~^M*, the median value indicated by observations of solar-type classical T Tauri stars as 



well as young very low- mass stars and brown dwarfs ( Andrews &: Williams|2005 2007 Scholz 



et al. 2006; Mohanty et al. 2012) 



Our choices of Rm, Rout and define the normalisation for any given surface density 
profile. Here we examine a scaled version of the Minimum Mass Solar Nebula (hereafter 
MMSN; H81), with a radial surface density profile given by: 



Sir! 



with So = 7.1 X lO^M^ gmcm" 



(13) 



H81 gives a somewhat larger value of 1.7x10'^ gm cm~^ for the normalisation Eq (for the 
Sun), because he uses essentially the same disk to stellar mass ratio as we do but a smaller 
outer disk radius {Rout ~ 36 AU instead of 100 AU). Our volume density is then 



p[r, z\ 



Por-^^/^e- 



(14a) 
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where the scale-height Zh[r] is given by equation (11), and 

po = 4.7 X 10-^° M. I ^^^^ I gmcm-=^ (146) 

Note that our expressions for po and Sq include an extra factor of M* compared to H81's 
formulation, due to the explicit dependence of our disk mass on the stellar mass. 



4.2. Ionization Rate and Chemistry 



We consider 3 sources of disk ionization: X-rays from the central star, cosmic rays, and 
the decay of radionuclides in the disk. For the first, we estimate the total X-ray luminosity 
{Lx) from the empirical median relationship Lx/L^, ~ 10~^'^ observed all the way from 
classical T Tauri stars to accreting brown dwarfs ( [Giidel et al.||2007 and references therein), 
and a suitable choice of L^, (described in §5). We consider coronal X-ray emission from an 
optically thin plasma dominated by impact excitation and ionisation by thermal electrons. 
Contributions from bound-bound, bound-free, and free-free radiation are included in the 
computation of our synthetic spectrum for all elements with atomic numbers 1 to 30, assum- 



ing the solar elemental abundances of Grevesse & Sauval (1998). We use line and continuum 



emissivities from the CHIANTI compilation of atomic data (Landi & Phillips 2006 and ref- 



erences therein), together with ion populations from Mazzotta et al. (1998), as implemented 
in the PINTofALE IDL software suite3 ( Kashyap fc Drake||2000 ) . The spectra are calculated 
using an isothermal X-ray emitting gas of temperature log(Tx) = 7.2, spherically distributed 
at a radius of 2i?^, 

X-ray ionization rates in the disc are subsequently computed with the moccasin code 



(Ercolano et al. 2003, 2005, 2008b), a 3-D photoionisation and dust radiative transfer code 



with a Monte Carlo approach to the propagation of ionising and non-ionising radiation. 
MOCCASIN has already been applied successfully to the modeling of X-ray irradiated proto- 



planetary discs (Ercolano et al. 2008a, 2009 Owen et al. 2010); more details regarding the 



physical processes and atomic data adopted for the calculations are given in Ercolano et al. 



(2003, 2008b) and Ercolano & Storey (2006). The first computations to address these issues 



were made by IG99, who also employed Monte Carlo methods to investigate the transport 
of X-rays within an idealised disk. Here we have repeated their calculations taking advan- 
tage of our newer code, the greater computing power available today, as well as a better 



^This ignores the harder X-ray component observed in some cases, e.g., dming flares; we intend to explore 
the additional effects of this hotter component in subsequent papers in preparation. 
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understanding of the X-ray properties of YSOs granted by results from the Chandra and 



XMM- Newton space telescopes available over the last decade (Getman et al. 2005 Giidel 



et al. 2007). Moreover, mocassin allows the inclusion of dust grains, the redistribution of 



secondary electrons and a full treatment of Compton scattering, and also uses a more so- 
phisticated Monte Carlo estimator for the radiation field. A detailed comparison of our rates 
to those of IG99 is the subject of a forthcoming paper (Ercolano et al 2012), where new pa- 
rameterised curves for use in, e.g., MRI calculations, are also provided. Here we simply note 
that for the same (grain depleted) elemental abundances, the ionisation curves we find and 
those computed by IG99 share the same physical characteristics, wherein a high ionisation 
rate is maintained in the upper tenuous disk regions before sharply falling off as the optical 
depth for absorption becomes of order unity, but the rate of decrease in our calculations is 
lessened by increased Compton scatterings due to the increased density. 



Additionally, we follow Umebayashi & Nakano (2009, hereafter UN09) in calculating 
the cosmic ray ionization rate for molecular hydrogen as: 



Hi 

CRJSM 



Sa[r,z] 



1 + 



Sa[r,z] 



3/4- 



-4/3 



+ e 



^CR 



Sb[r,z] 



^CR 



3/4" 



-4/3- 



(15) 

where Sa[r, z] = p[r, z] dz is the surface density of the disk above the height z at radius 
r, and ^^[r, z] = p[r, z] dz is similarly the surface density below the height z. Cch ism — 
10^^'' s~^ is the interstellar cosmic ray ionization rate for molecular hydrogen, and T^cr ~ 
96 g cm~^ is the characteristic attenuation surface density for cosmic rays (UN09). The two 
right-hand side terms in equation (15) thus express the total cosmic ray ionisation rate at a 
location within the disk as the sum of the rates due to attenuated cosmic rays arriving from 
the hemisphere above and the hemisphere below that location. 



Finally, we also follow UN09 in calculating the ionisation rate of molecular hydrogen 
due to the combined effect of short- and long-lived radionuclides as 



7.6 X 10-^^ s-\ 



The total ionization rate ( of H2 is the sum of the X-ray, cosmic ray and radionuclide 
rates: C = + Ccr + Cr^ ■ Fo^' the ionization balance, we use a simplified chemical network, 
including both gas-phase and grain chemistry, essentially described by model. 4 of 1N96. 



The gas-phase chemistry in this model is that adopted by Oppenheimer & Dalgarno (1974 
hereafter OD74), and schematically described by: 
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Table 1: Gas Phase Chemistry (Reactions and corresponding Rate Coefficients) 

+ e~ ( (ionization rate, s~^) 

a = 3 X 10~VVT cm\-i 
+ hu 7 = 3 X W^^/VT cm^s"^ 
+ M+ ;5 = 3 X 10"^ cm^s-i 
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where m and m+ represent a molecule and its ionized counterpart, M and M+ are a neutral 
heavy metal atom and its ionized counterpart, and e^ is a free electron in the gas phase. 
is the total ionization rate for hydrogen molecules, as calculated above; the temperature for 
the reaction rates a and 7 is given by our radial temperature profile, equation (9). 

In our case (as in IN96), the metal is magnesium and the neutral molecule m is H2. 
However, the important ionized molecule m"*" is not ; while the ionization of H2 initially 
creates this ion, the subsequent reaction chains invariably result in the dominant m+ being 
an anion containing both hydrogen and some other heavier element (s), such as HCO+, NH^+ 
etc. Following OD74, we only consider species involving hydrogen, carbon and oxygen for 
the molecular reactions, and omit nitrogen; we further omit those OD74 reactions involving 
pure atomic or molecular oxygen, following PCll (see latter for justification). In this case, 
the reactions Rl and R2 in Table 1 may be explicitly expanded into the following chain: 

H2 + XCRR ^ H+ + e- R(a) 

H+ + H2 ^ H+ + H R{b) 

H+ + CO ^ HCO+ + H2 R(c) 

HCO+ + e- ^ CO + H R(d) 

2H + gr ^ H2 + gr R(e) 

where XCRR signifies either an ionizing X-ray, cosmic ray or particle released by radionuclide 
decay, and gr represents a grain surface which catalyses the conversion of atomic to molecular 
hydrogen. The cascade R(a)-(c) boils down, without any loss of generality, to: 

H2 + XCRR + CO ^ HCO+ + e- + H R(/) 

which shows that the dominant molecular ion here is HCO'^. Reactions R(/), (d) and (e) 
capture the essence of our gas-phase chemistry: in the presence of ionizing infiuences and 
CO, molecular hydrogen is converted to HCO"*", electrons and atomic hydrogen; the HCO+ 
and electrons dissociatively recombine to give back CO and more atomic hydrogen; and the 
circle is closed by the hydrogen atoms coalescing into hydrogen molecules on on the surface 
of grains that are well-mixed with the gas. 
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However, this set presents a problem when the abundance of well-mixed grains becomes 
extremely low. In such circumstances, the rate of conversion of H to H2, reaction R(e), can 
decrease to the point where atomic hydrogen dominates over molecular in the disk (which 



is indeed what happens in the dust-free simulations by Bai & Goodman (2009, hereafter 
BG09): X. Bai, pvt. comm. 2011; see detailed discussion in §6.1). What is the threshold grain 
abundance below which this occurs? We assume grains of uniform radius 0.1 /xm (see end of 
this section), and adopt the reaction rate formula supplied by BG09 (their equation (27)). 
Then, for a fiducial disk temperature of 100 K (achieved at ~15 AU around a solar-type star 
and 5AU around a very low-mass one), a fairly low gas number density of uh ~ 10^ cm^'^ 
(attained at 2-3 scale heights at these radii), a sticking coefficient of Sh ~ 1, and a constant 
H-to-H2 conversion efficiency of e ~ 10~^ (value preferred by BG09 after correcting for an 
error in the formulae by Cazaux & Tielens (2004)), the timescale for converting initially 
completely atomic H entirely into H2 reaches ~10^yr (the maximum lifetime of disks), for 
a dust-to-gas mass ratio of Rd:g = 10~^. Thus, we nominally expect atomic hydrogen to 
dominate over H2 in the upper/outer (low density) parts of a disk, when R^ g < 10~^. 

There are however a number of uncertainties in this calculation, most importantly in 
the conversion efficiency e. On the one hand, Cazaux & Tielens (2002, 2004) show that e 
should increase rapidly with decreasing temperature (i.e., in the outer disk), compared to 
the fixed value adopted by us and BG09, thereby lifting the H2 fraction. On the other hand, 
we (in common with previous works such as BG09, Bll and PC 11) have not considered the 
formation of ice mantles on grains beyond the ice line (i.e., at T^isk ^ 170 K). H cannot be 
chemisorbed onto icy surfaces, only physisorbed; since such physisorption ceases at >20 K, 
H2 formation on icy grain mantles halts at these temperatures (see discussion in Cazaux & 
Tielens 2002). Conversely, H2 formation under such conditions may be mediated by other 
molecular species (Cazaux & Tielens 2002, and references therein). 

Given these uncertainties, we assume that the gas always exists mainly as H2, in accor- 
dance with the standard view of protoplanetary disks. For the two higher dust abundance 
cases we examine, Rd-g = 10^^ and 10"^ (see end of this section), this assumption is con- 
sistent with the H-to-H2 conversion rate discussed above. For our two very low abundance 
cases, Rd:g = 10~^ and 0, this amounts to assuming the presence of channels other than 
grain-surface recombinations for the conversion. Imposing H2 dominance by fiat, for all dust 
abundances, allows us to focus here on the primary process by which dust depletion affects 
the MRI: as grains are removed, the rate of electron adsorption onto grains (and the at- 
tendant rate of ion-electron recombinations on grains) nosedives, strongly altering the MRI 
efficiency. If our calculations indicate that very low grain abundances are physically relevant 
(as turns out to be the case when we examine accretion rates; see point (v) in §6.3), it will 
serve as an impetus to investigate the additional effect of hydrogen recombination channels 
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at such abundances in more detail in future work. We note that Bll similarly impose H2 
dominance in their grain-free case, for the same reason (discussed further in §6.1). 

To enforce H2 dominance, we approximate the gas-phase reactions R((i)-(/) by the 
following set in all cases (both with and without grains): 

H2 + XCRR HCO+ + e" R(^) 

HCO+ + e" ^ H2 R{h) 

We also use ^{h) to approximate the dissociative recombination of HCO"'' and e~ on grains, 
if the latter are present (we do account for the difference in recombination rates between 
the gas phase and grain surfaces). Despite first appearances, this system preserves exactly 
the number of hydrogen particles: both the formation and dissociative recombination of an 
HCO"'" ion liberate an H each, which are simply lumped together as one hydrogen molecule 
in R(/i), expressing our condition that H2 be the final state of any free hydrogen. Reactions 
R((yf)-(/i) allow us to accurately follow the equilibrium abundance of the dominant molecular 
ion HCO^, while enforcing the conversion of H into H2 regardless of the grain abundance. 
We do not follow the destruction and reformation of CO in this network, since HCO+ is a 
trace species relative to CO. 

In summary, the gas-phase reactions Rl and R2 are interpreted as R(5') and (/i), with 
the stipulation that m = H2 and m+ = HCO^; their rates remain as given in Table 1. 

As clear from the above, the chemical model further includes dust grains, and chemistry 
due to gas-grain as well as grain-grain collisions. Grains are assumed to be either neutral, 
singly charged or doubly charged (henceforth grO, gr^*" and gr^^*" respectively)]^ Reactions 
involving grains are divided into two classes: mantle chemistry, in which a neutral gas- 
phase particle is adsorbed onto a grain mantle (or the reverse process, desorption, in which 
an adsorbed neutral species is released into the gas phase); and grain chemistry, wherein 
a charged particle (whether gas-phase or grain) collides with a grain resulting in charge 
transfer. We refer the reader to IN06 for the specific mantle and grain reactions considered 
(their Tables 3 and 4) and an exhaustive discussion of the model (reaction rates, sticking 
probabilities, desorption etc.). 

Finally, the set of equations is closed by specifying the total number density of particles 
in the gas-phase, n; the abundance of the heavy metal species (Mg, both neutral and charged) 



^Ideally, higher grain charges should also be considered. In reality, however, this is unlikely to make a 
significant difference to our results for 0.1 sized grains: our derived electron abundances are in excellent 
agreement with those of PCll for the same grain size (see discussion in §6.1, point (vi)), even though the 
latter authors allow a much larger range of grain charges. 
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relative to hydrogen, x[Mg]; the density and radius of each grain; and the dust-to-gas ratio. 
The chemical reactions cause negligible variations in the mean molecular weight fi of the gas 
(IN06), so we adopt a constant yU = 2.34, corresponding to a solar composition for the most 
abundant elements, H and He. The number density of gas-phase particles at any location is 
then n[r, z] ~ p[r, z]/ {fiuiH), with p given by equation (14) and mn the mass of a hydrogen 
nucleus. For the gas-phase abundance of Mg we adopt x[Mg] = 6.78 x 10~^ per H2 particle 
(10-^ solar). 

Real disks have a distribution of grain sizes, determined by the balance between colli- 
sional coagulation and fragmentation. Recent simulations by Birnstiel et al. (2010) indicate 
that fragmentation results in a significant population of small grains even in the face of 
strong grain growth. Moroever, while larger grains quickly settle to the disk midplane, very 
little turbulence is needed to keep the smallest ones aloft and well-mixed with the gas. Over- 
all, therefore, one expects the bulk of the dust mass, locked in larger grains, to be rapidly 
confined close to the midplane, while a significant fraction of the grains in terms of number, 
skewed towards the smallest sizes, remains mixed throughout the disk. Finally, the smallest 
grains are the most important in the grain reactions discussed above: the latter depend 
crucially on the grain surface area per unit gas mass, which is largest for the smallest grains. 

We encapsulate this picture as simply as possible in our calculations, as follows. We 
consider only very small ISM-sized grains, of uniform radius a = 0.1 fim and internal density 
Pgr = 3 g cm~^, well-mixed with the gas. At the same time, we vary the dust-to-gas mass 
ratio from run to run, using the values 1:100 (standard ISM), 1:10'^ (90% depleted relative 
to standard ISM), 1:10^ (99.9% depleted) and (100% depleted). The "depleted" cases 
simulate conditions wherein the dust mass is increasingly bound up in larger grains that 
have settled to the midplane and thus no longer interact with the gas. The 100% depletion 
case is not expected to hold literally in real disks, but rather represents the limiting value. A 
more realistic treatment demands dynamic calculations that account for the time-dependent 



effects of grain growth and settling (Turner et al. 2010); this will the subject of future work 



The methods we use to find chemical equilibrium are entirely standard. We time- 
integrate the chemical kinetic equations by semi-explicit extrapolation, an algorithm that 



works well with such stiff ODEs and is described in Numerical Recipes (Press et al. 1992). 
The integration is stopped when all species' timescales for further change become much longer 
than the disk lifetime. Substituting the final abundances back into the kinetic equations in 
a few specific cases demonstrates that the solutions are very close indeed to equilibrium. 



4.3. Resistivities 



Armed with the equihbrium abundances of the various species from the chemical network 
calculations, we compute the resistivity at every location, and thereby investigate where the 
disk is MRI active by the criteria of §2 for a given field strength B. We follow Wardle (2007) 
in writing the Ohmic, Hall and Pedersen conductivities (ao, cth and ap respectively) as 

j 

ec sr^ tlaZ. 



ec V ^ nj\Zj\l3. 

ap — 



2 



where the summation is over all charged species j with particle mass m^, number density rij 
and charge ZjC. The Hall parameter /3j (not to be confused with the plasma (5 parameter) is 
the ratio of the gyrofrequency of a particle of species j to its collision frequency with neutrals 
(of mean particle mass m„ = ^itlh and density p„): 

\ZAeB 1 , , 

13, = ^ 19 

Here jj =< av >j /{rrij + m„) is the drag coefficient and < av >j the rate coefficient for 
coUisional momentum transfer between species j and the neutrals, making 7jp„ the collision 
frequency with neutrals. This allows the induction equation to be written in the form: 

_ = Vx(vxB)-Vx [7]o{V X B) + r]H{V x B) x B + r]A{V x B)^] (20) 

where denotes a unit vector and " _L " indicates the component of a vector perpendicular 
to B. The first expression on the right is the inductive term (J), while the second, third 
and fourth are the Ohmic (O), Hall (H) and ambipolar (A) resistive terms respectively. The 
ratio of the inductive term to whichever of the resistive terms is dominant is thus of the form 
Lv/f] ~ L'^/tt], where L is the field diffusion lengthscale and t the field advection (growth) 
timescale, and the resistivity rj equals rjo, rju or 77^; for the field to grow, one expects this 
ratio to exceed unity. L ~ vk/^ for large-scale fields generated by shearing radial ones, 
while L ~ Vji^/Q for local fields created by the MRI; in both cases t ~ This yields 

equations (1) and (2) when Ohmic diffusion holds sway {1] = rjo), and the naive condition 
equation (3) when ambipolar diffusion rules instead [rj = rjA)- The resistivities are: 
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.« = T^^ (22) 

47ro"_L cr^ 

= Vo (23) 

47rO'_L 0"_L 

with (Tx = a/o"!^ + (Tp being the total conductivity perpendicular to the magnetic field. Note, 
from equations (16)-(19) and (21)-(23), that: (i) rjo is independent of the field strength B, 
while riH and rjA scale linearly and quadratically, respectively, with B; and (ii) if ions and 
electrons are the the only charged species {not true in our case), then equation [25] reduces 
to rjA = f^if^eTjo, yielding (using /3j -C f3e) the simplification v^z/va^ ~ liPi/^, as advertised 
in §1.2. To compute the resistivities, we must specify the rate coefficients < av >j for j = e 



(electrons), i (ions) and gr (grains). We employ the values supphed by Wardle & Ng (1999): 



128A;r, 



1/2 



<av>e= 10"^^ crn^ s"^ (24) 

V dnnie J 

< av >i = 1.6 X 10"^ cm^ s"^ (25) 

2 / \28kTn\ ^ 3—1 



< av >„r = TTa^ cm^ s"^ (26) 

where Te and T„ are the electron and neutral temperatures respectively, and a is the grain 
radius (fixed at 0.1 /xm in our case). Note that < av >i is the same for molecular and metal 
ions. We assume that Te = T„ = T[r], with the latter given by equation (9). 



5. Stellar Models Considered 

We examine two stellar models. The first corresponds to a typical classical T Tauri star: 
M^, = O.TMq , R^: = 2Rq and T^, = 4000K. The second corresponds to a very low mass star, 
also in its disk accretion phase: M^, = O.IM© , R^, = IRq and T, = 3000K. The radii and 



temperatures adopted for both masses are from the Lyon evolutionary tracks (Baraffe et al. 



1998) for an age of 1 Myr, appropriate for objects in their classical T Tauri phase. 



These numbers correspond to bolometric luminosities of L^, = 0.92 and 0.07 Lq for M^, 
= 0.7 and 0.1 Mq respectively. Our adopted scaling Lx/L^ = lO^'^'^ then yields an X-ray 
luminosity of Lx = 1.12 x 10=^° erg s"^ for O.7M0, and Lx = 8.84 x lO^^ erg s'^ for O.IM©. 
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6. Results 

6.1. Ionization Balance and Chemistry 

The equilibrium fractional abundances of all the species we have considered, defined 
as Xa = Tia/nn^ (where and riHa are the number densities of species a and H2 molecules 
respectively), are plotted in Figs. 1 and 2. Fig. 1 shows the situation at a radial distance of 
1 AU, for stellar masses of 0.7 Mq (top panels) and 0.1 Mq (bottom panels), and dust-to-gas 
ratios (assuming 0.1 /im grains) ranging from 10~^ (standard ISM) to zero (left to right). 
Fig. 2 shows the same for 10 AU. Some clear trends are noteworthy, as follows. 

(i) At any radius, the inclusion of even a tiny fraction of well-mixed 0.1 /im grains drastically 
suppresses the abundance of free electrons as one moves towards the midplane, relative to 
the no-dust case. Moreover, the height at which the drop-off commences and the steepness 
of the decline both increase with grain abundance. In other words, in the presence of grains, 
Xe decreases dramatically with both increasing density and increasing grain abundance. 

(a) In low density regions of the disk, i.e., sufficiently high above the midplane, charge 
equilibrium is dominated by the balance between free electrons and molecular ions (HCO"'" 
in our case), independent of the grain abundance. 

(Hi) The abundance of metal ions, XMg+, falls off precipitously towards the midplane with 
increasing grain abundance. 

(iv) Related to the above, we see that in the absence of grains, charge equilibrium close to 
the midplane (i.e., at high densities) is maintained by the balance between metal ions and 
free electrons. At sufficiently high grain abundances, on the other hand, metal ions become 
negligible, and charge equilibrium near the midplane is maintained instead by either the 
balance between charged grains, or between free electrons and molecular ions (HCO"'"). 



BG09 and Okuzumi (2009) have previously noted the above trends for solar-type stars; 
our analysis reproduces their results for such stars {^^0.7 Mq), and shows that these trends 
apply to disks around very low-mass stars (~0.1 Mq) as well. 

(v) The reason for the decrease in metal ions with increasing grain abundance has been 
debated before. IN06 argue that the decline is due to the sweeping up (via adsorption) of 
Mg atoms by grains, preventing the formation of Mg+. In contrast, BG09 posit that it is 
because grains promote the recombination of metal ions into neutral atoms. Our analysis 
shows that both effects can be at play, depending on the disk conditions, as follows. 

To begin with, note that at 1 AU for 0.7 Mq, the abundance of free (gas-phase) neutral 
Mg atoms close to the midplane is hardly affected by the grain abundance, while the abun- 
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dance of Mg adsorbed onto grains (i.e., grMg) is orders of magnitude below that of the free 
neutral Mg. Clearly, most of the neutral Mg is not swept up by grains, and thus such adsorp- 
tion does not cause the steep decline in Mg"''. However, the abundance of adsorbed Mg does 
increase as the Mg"*" and electron abundances decrease. The only explanation is that grains 
facilitate the recombination of metal ions into neutral atoms, which are then adsorbed onto 
the grains. In other words, the cause of the decline in Mg+ here is recombination stimulated 
by grains, and the increase in adsorbed Mg is a result of this process. This is consistent with 
the scenario proposed by BG09. The same is true at 1 AU for 0.1 Mq as well (the abundance 
of adsorbed Mg is much higher than in the 0.7 Mq case, but still far below that of neutral 
Mg, which remains unaffected near the midplane by the grain abundance). 

At 10 AU, the situation is radically different. As the grain abundance increases, the 
fractional abundance of adsorbed Mg becomes comparable to, and eventually much higher 
than, the abundance of gas-phase neutral Mg (at dust-to-gas ratios of > 10~^ for 0.7 M©, 
and > 10~^ for 0.1 M©). At the lower temperatures prevalent at this radius, grains do sweep 
up much of the neutral Mg, directly preventing the formation of Mg"*", as IN06 propose. This 
effect is also much stronger for the 0.1 Mq case compared to 0.7 Mq. 

(vi) Our quantitative electron abundances in the grain-free case agree with the results of 
PCll and BG09, but differ significantly from those of Bll, for a subtle reason. 

To illustrate this, we first plot and the ionization rate C as a function of the total 
number density of particles, at lAU for a 0.7 M© star, for both the grain-free case and a 
standard dust-to-gas ratio of 10~^ (Fig. 3a). For comparison, we refer to Fig. 12 by PCll. 

In the latter, PCll plot function of the number density of H2 (essentially equal to 

the total number density), for a fixed disk temperature of 280 K (corresponding to 1 AU for 
a solar-type star in the H81 model) and a constant C, = 10^^^ s^^. They show the results for 
both their own chemical network as well as the simple and complex networks of BG09. 

Our Fig. 3a shows that ( = 10^^^ s^^ corresponds to n 2 x 10^^ cm^^ in our model at 
1 AU. For this density, we find Xe ~ 10~^^ for a gas-to-dust ratio of 10~^ with 0.1 /im grains, 
and ~ 3 X 10~^° without grains. This is in excellent agreement with the x^. plotted by 
PCll for the same density and same grain conditions. Specifically, our Xg almost exactly 
equals the PCll network result at this density for both the dust- free and dusty cases; it also 
equals the BG09 results for the dust-free case, and is within a factor of ~3 of the latter for 
the standard dust-to-gas ratio. The agreement with PCll is unsurprising, given that our 
chemical networks are nearly identical. Importantly, this also shows the negligible impact of 
the difference in disk temperature between our models (393 K versus 280 K at 1 AU). 

Next, we compare to the results of Bll. To make this straightforward, we replot our Xe 
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and C at 1 AU for the 0.7 Mq case as a function of fractional scale height {z/zh) in Fig. 3b, 
and compare to the analogous plot in Fig. 1 by Bll. Otii ionization rates due to cosmic rays 
and radionuclides differ somewhat from Bll's, so it is more instructive to compare our results 
in the regions where X-ray ionization dominates, i.e., > 10^^^ s^^. For specificity, consider 
again the exact value = 10^^^ s^^, reached in both our case and Bll's at z/zh ~ 2 (marked 
in our plot as well as Bll's by a vertical dashed line). For the case of 0.1 /xm grains with a 
dust-to-gas mass ratio of 10~^, Bll finds 10 at this (, roughly comparable to 

our value of 10~^^. For the dust-free case, however, he finds Xe ~ 10~^^, more than 2 orders 
of magnitude below our value of 3 x 10~^°, and thus also similarly lower than the PCll and 
BG09 predictions. The difference appears to be due to a subtle change in chemistry in the 
grain-free case between the earlier analysis by BG09 and the new one by Bll (the following 
discussion was greatly aided by X. Bai, private comm., 2011). 

In particular, BG09 and Bll differ in their treatment of H2 formation. Some of the 
intermediate chemical reactions lead to the formation of atomic hydrogen, requiring a channel 
to convert this back to H2 (see discussion in §4.2). BG09 assume that the recombination of 
atomic to molecular hydrogen occurs on grain surfaces. In the absence of grains, however, 
most of their hydrogen remains in atomic form. Bll, aware of this issue, artificially adds 
a reaction to convert H to H2 independent of grains (X. Bai, pvt. comm. 2011), to ensure 
that H2 is the dominant form in all cases. This in turn increases the NH4 abundance in his 
chemical network by orders of magnitude compared to BG09, and thereby decreases the free 
electron fraction by a comparable amount due to the dissociative recombination of NH4 . 
This explains the difference between BG09 and Bll. 

Conversely, our calculations (and those of PCll) have only HCO"*" and no NHj , but 
we also enforce by fiat (essentially as Bll does) that any atomic hydrogen is converted to 
H2 even for the grain- free case (see §4.2). We therefore see that the combination of mainly 
atomic hydrogen and the inclusion of NH4 serendipitously leads BG09 to approximately 
the same Xe as we (and PCll) get with mainly molecular hydrogen in the absence of NH^, 
but the combination of molecular hydrogen with NH^ leads Bll to much lower Xg. A full 
resolution of this issue requires a better understanding of hydrogen recombination channels 
for extremely low grain abundances. 

Nevertheless, as noted earlier, even a small abundance of well-mixed grains strongly 
affects the electron fraction. In the presence of dust, as shown, our results are consistent 
with those of Bll as well as BG09 and PCll, despite variations in the chemistry. 
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6.2. MRI Active Zones 

Next, we calculate the location and vertical extent of the MRI-active zone at every 
radius, using the stress-maximizing technique described in §3. Before discussing our full disk 
results, we illustrate the technique for the fiducial case of 0.7 Mq with a dust-to-gas ratio of 
10~^, at disk radii of 1 AU and 10 AU (Figs. 4a and b respectively). 

6.2.1. Illustrative Results at 1 and 10 AU for 0.7 Mq and Standard Dust Abundance 

Consider the 1 AU case first. The top left panel of Fig. 4a shows hB"^^ where h is 
the thickness of the active layer, as a function of the field strength B. We see that this 
quantity increases monotonically with and then drops extremely rapidly to zero beyond 
its maximum value, as 13 falls below ^min for stronger fields and the MRI is quenched at all 
vertical locations at this radius. In other words, the maximum value of hB"^ occurs essentially 
at the Smax allowed for a given radius; hence, for all practical purposes, we are justified in 
setting Smax (which corresponds to I3min) precisely equal to the (very slightly lower) field 
strength at which the height-integrated stress hB"^ is maximized, as stated in §3. 

The top right panel of Fig. 4a explicitly shows the vertical location of the MRI-active 
region as a function of the field strength, and the conditions that constrain its thickness. The 
horizontal dashed line marks the particular value B^^, at which hB"^ is maximized in the 
previous panel. We see that the upper envelope of the active layer is always determined by 
the limiting ambipolar condition (5 = /Smm- Conversely, for small B, the lower envelope is set 
by the minimum Ohmic Elsasser number criterion, A_4 = 1. As the field strength increases to 
nearly -Bmax, however, the lower envelope of the active region instead becomes determined by 
the ambipolar criterion (3 = /3min as well. This is not surprising: given that for fields stronger 
than Smax, the ambipolar condition shuts off MRI everywhere (by the definition of Smax), it 
is to be expected that, close to -Bmax, this criterion supersedes the Ohmic one even close to 
the midplane. This is ultimately because the ambipolar resistivity tja increases quadratically 
with B, while rjo is independent of B. While Bll does not explicitly mention this effect, it 
is present in his grain-inclusive calculations as well, as evinced by the characteristic "dome- 
shaped" peak of the active region close to -Bmax in both his plots and ours (see his Fig. 3, 
for 0.1 yum grains at 1 and 10 AU), arising from the ambipolar condition sculpting both the 
upper and lower envelopes of this region. Finally, the vertical line shows the height at which 
the magnetic Reynolds number criterion for large-scale fields achieves its limiting value, 
Ak — 10. Note that this location is independent of the field strength. Moreover, for any 
value of B for which an active layer exists, we see that the Elsasser number criterion is 
achieved higher in the disk (i.e., at smaller rjo) than the Reynolds condition, which refiects 
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the fact that vk ^ for a thin disk in which ^ < /^min < 1 (see discussion in §2). 

The bottom left panel of Fig. 4a shows the values of the MRI-activity criteria, Kk, 
and P/^min, as a function of height for the optimal field strength -Bmax- The active 
region (A^ > 10, A^ > 1 and ^/^mm. > 1) is sandwiched between z/zh ~ 2.5 and 3.5 (the 
intersection of the B — Bj^ax line and the active zone in the previous panel). As noted 
above, we also explicitly see that both the upper and lower limits of this layer are set by the 
ambipolar criterion: the Ohmic Elsasscr number exceeds unity slightly deeper in the disk 
than the actual bottom of the active region. Finally, the far less stringent Reynolds criterion 
is satisfied much deeper in the disk. Regions closer to the midplanc than the locations of 
both the Reynolds and Elsasscr limiting values comprise the Dead Zone, regions that lie 
between the two make up the Undead Zone, and regions where the ambipolar criterion is 
not satisfied (everywhere outside the active layer in this case) constitute the Zombie Zone. 

The bottom right panel of Fig. 4a shows the Ohmic, Hall and ambipolar resistivities as 
a function of height, for -Bmax- We see that the curves for A_4 and /3//3min, in the previous 
panel, essentially reflect the inverse of r)o and tja, as expected, rjo, which depends on the 
ionization fraction, decreases with increasing height (as the electron abundance increases), 
while TjA, which depends on both the ionization fraction and the density, first decreases with 
increasing height (as Xe increases) and then rises again (as the density continues to decrease). 
Note that rjH changes sign at z/zh ~ 5. Note also the flattening of all the resistivities below 
z/zh ~ 2, which is caused by a flattening of the ionization rate close to the midplane (where 
the ionization is mainly due to cosmic rays and radionuclides, at this radius; e.g., see Fig. 3b). 
The slight jitter in t/a very close to the midplane is due to small number photon statistics, 
at these very high optical depths, in our Monte Carlo transfer simulations. 

The results for 10 AU, plotted in Fig. 4b, are qualitatively very similar, albeit with three 
differences worth mentioning. First, the range of flcld strengths for which an MRl-active 
layer exists, and thus -Bmax, is signiflcantly smaller than at lAU. This is because the gas 
pressure Pgas at all heights at this larger radius is lower, and thus the magnetic pressure 
Pb for which /3 > ^min is also smaller. Second, Ak > 10 at all heights at 10 AU (which is 
why the curve for this criterion does not appear in the second and third panels of Fig. 4b), 
meaning that there is no Dead Zone at this radius, only Undead and Zombie Zones. Third, 
there is no flattening of the 7] curves close to the midplane, unlike in the 1 AU case: the 
lower densities at 10 AU allow X-rays to penetrate much more easily to the midplane, and 
thus dominate over cosmic rays and radionuclides in controlling the ionization rate there. 

These results are very similar to those of Bll for the same stellar type, disk radii and 
grain properties. The only small difference is that our -Bmax arc somewhat greater than his (by 
3-5 x) at a given radius, because of differences in chemistry, and also because our ionization 
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rates are somewhat higher for the same ionizing flux (due to a shghtly lower surface density 
in our case compared to the standard MMSN - sec discussion of equation (13) in §4.1 - and 
also due to our inclusion of Compton scattering in a self-consistent manner, which makes our 
X-ray ionization rate fall off slower with increasing optical depth compared to the results 
of IG99, which Bll uses for his analysis). Conversely, as discussed in §6.1, our electron 
abundances close to the midplane for a dust-free case are considerably higher than Bll finds; 
for this specific idealized situation, our B^^x are correspondingly significantly higher (not 
plotted), because both rjo and tja are smaller when more free electrons are present. Finally, 
our results for the 0.1 Mq case are qualitatively similar to, with the same trends as, the ones 
plotted in Fig. 4 for 0.7 Mq. 



6.2.2. Full Disk Results 

With these preliminaries dealt with, we move on to a discussion of our full disk results. 
In Figs. 5-8 we plot our results for M* = 0.7 Mq, for dust-to-gas ratios of 10~^-0; Figs. 9-12 
show the same for 0.1 Mq. The field strength adopted at each radius is the optimal -Bmax 
derived as discussed above. The top panel in each plot illustrates the relative importance of 
the various resistivities as a function of radius and column depth, the middle panel addition- 
ally shows the locations of the Dead and Undead Zones, and the bottom panel shows the 
location of the Zombie Zone as well. The final MRI-active layer in each case is the region 
bounded by black filled circles in the bottom panel. We note that there are a few sharp dis- 
continuities/jitters in the boundary of the active layer at large radii in some cases, marking 
a transition in the shape of the layer with radius. The changes in the shape are themselves 
real; they reflect the fact that the resistivities vary with radius and height, sometimes very 
rapidly (as illustrated in Figs. 1 and 2). The few discontinuities and sharp fluctuations at 
these transitions, however, are not real, but rather due to the inherent statistical uncer- 
tainty in our Monte Carlo determination of the ionization rate. Fewer photon packets are 
incident on the disk at large radii, yielding a relatively higher jitter in the derived ioniza- 
tion rates there. The resulting small statistical oscillations in the inferred resistivities can 
lead to somewhat larger fluctuations in the active layer boundary, by pushing the diffusivity 
over the critical threshold that deflnes the boundary. We have made no attempt to remove 
these jitters, since they are a measure of the accuracy of our results (and smoothing them 
out by using more X-ray packets would have entailed prohibitively long processing times). 
Nevertheless, the overall shape of the active layer as a function of radius is clear in all cases, 
and the similarity in shapes between the 0.7 and 0.1 Mq simulations for a fixed dust-to-gas 
ratio (and the agreement with the semi-analytic results of Bll for the dusty solar-type case, 
discussed above) implies that our results are physical, and not dictated by statistical errors. 
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There are five main conclusions that we may draw from these plots. 

(i) For both 0.7 and 0.1 Mq, the Dead and Undead Zones become smaller with decreasing 
grain abundance, and essentially disappear in the dust-free case. This is a consequence of 
increasing Xg, and thus decreasing Ohmic resistivity, with declining amounts of dust. 

(a) For a fixed dust-to-gas ratio, the radial extents of the Dead and Undead Zones are smaller 
for 0.1 Mq than for 0.7 Mq, by a modest factor of ~ 1.5-2. The X-ray ionizing flux from the 
very low-mass star is far lower than that from the solar-type one, but the gas density (and 
hence X-ray optical depth) at any radius is also much lower in the disk around the former; 
the two effects compensate for each other in setting the ionization rate, making the Dead and 
Undead Zones roughly comparable in size between the two stars. For the solar-type case, 
our Dead and Undead Zone extents (e.g., ~3 and 13 AU respectively for a dust-to-gas ratio 



of 10 ) are also consistent with those found by Turner & Drake (2009, who use somewhat 



larger 1 /zm grains, and a magnetic pressure fixed at 0.1% of the midplane gas pressure). 

(in) In all cases, both the upper and lower boundaries of the MRI-active layer are set by the 
minimum ambipolar condition /3 = /3min- At radii where the Undead Zone exists, its upper 
envelope (set by the Elsasser number limit A^ = 1) roughly follows, but still lies below, the 
lower envelope of the active layer. This has an important consequence for mass accretion. 
The raison d'etre for the Undead Zone is that the MRI is not required for angular momentum 
transport within this zone: the transport here occurs via large-scale toroidal fields generated 
from seed radial ones by Keplerian shear, as long as the Undead Zone is able to siphon off 
seed fields from neighboring MRI-active areas. If this zone is separated from the MRI-active 
region by an ambipolar-dominated inactive layer, however, such siphoning may be quenched. 
In this case, the Undead Zone effectively dies as well, and cannot be invoked to supplement 
the MRI-driven angular momentum transport and mass accretion. 

(iv) For a fixed dust-to-gas ratio, the shape and location of the active layer is roughly 
similar for 0.1 and 0.7 Mq. The midplane in the outer disk is MRI active for all the dust 
abundances considered, while activity extends to the midplane in the inner disk only in the 
dust-free case. For both stars, the active layer moves closer to the midplane with decreasing 
dust abundance (because the fraction of free electrons increases at all heights, and thus r]A 
decreases); simultaneously, it also becomes thinner, in units of decades of column density. 

(v) The combined effect of the two trends in (iv) is that the total mass of the active layer 
remains roughly constant with varying dust-to-gas ratio, comprising ~5-10% of the total 
disk mass for both 0.1 and 0.7 Mq. 
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6.3. Accretion Rates 

Having determined the locus of the active layer in each disk, we estimate the mass ac- 
cretion rate through it. As discussed in §3, Bll adopts M = hB'^^.^^/4Q, which only holds for 
a radially constant rate, while we use the general expression M = {l/2rQ) d{r^hBl^^^)/dr. 
Figs. 13-16 show our results for dust-to-gas ratios of 10~^-0. In each plot, the top panels 
illustrate various relevant quantities as a function of disk radius for 0.7 Mq, while the bottom 
panels show the same for 0.1 M©. For each star, we plot from left to right: (1) -Bmax; (2) 
^-^max) which determines the accretion rate (modulo fl) in Bll's case; (3) r'^hB'^^^, the slope 
of which determines the accretion rate (modulo rfl) in our case; and (4 ) the M derived from 
these quantities, using both Bll's formula and ours. As an aside, notice that the shape of 
the Bll M as a function of radius is similar to that of r'^hB^^^; this merely reflects the fact 
that the former is proportional to r^^^hB^^, only a factor of r^/^ different from the latter. 
We further note that we have evaluated the differential in our M expression via simple fi- 
nite differencing over our radial grid, which leads to small-scale unevenness in the derived 
accretion rate. We have not attempted to smooth this out, since the overall trends in our 
M - average magnitude as a function of radius, steep rise/falloff at certain locations, and 
a reversal in the direction of the flow at others, all of which are real, as discussed below - 
remain clear. There are six main implications of these plots. 

(i) The M inferred with Bll's formula are radially variable. While Bll discusses some phys- 
ical consequences of this, the bare fact of radial variations demonstrates the inconsistency 
in using this formula, which is only valid for M constant with r. 

(ii) Our M, derived from the general equation, is also radially variable; moreover, it also 
changes sign at some locations. That this must happen is easily verified by eye from the plot 
of r^hB^^^ in the third panel: where this quantity increases with radius, the accretion rate, 
proportional to the slope of this function, is positive; where it is fiat, the rate drops to zero; 
and where it decreases with radius, the rate is negative, i.e., material flows outwards. 

(in) The immediate consequence of these radial variations in M is that material builds up in 
some regions of the active layer, and is preferentially evacuated from others. Whether this 
process runs away, oscillates or reaches a steady state, however, is complicated, as follows. 

Consider, for example, a disk radius at which the active layer does not extend all the 
way to the midplane, and where moreover there is a net influx of matter (M coming into the 
region exceeds that going out). The latter enhances the density at every vertical location 
inside the active layer here. Therefore, at any given time, and for a fixed field strength, the 
new ionization and chemical equilibrium at each height within the layer here equals that at 



-27- 



a lower (denser) height at an earher tim^ But the bottom of the active layer is originally 
set by the density above which (or more precisely, the ionization fraction below which) the 
MRI cannot be sustained for 5max- Thus, with increasing density (and hence decreasing Xe), 
the bottom of the layer moves upwards, while the location of the top remains fixed for the 
given field strength (since there is no change in density at any point originally above the 
active layer). Nor can the field strength change from the initial -Bmax^ it cannot get weaker, 
because we have assumed that the MRI maintains the strongest possible field that still allows 
the MRI to operate, and it cannot get stronger, because -Bmax is already by definition the 
strongest (vertically constant) field that can be maintained here. The upshot is that the 
active layer keeps getting thinner, and thus the quantity hB^^^ continues to diminish, at 
this location; the absolute rate of change in hB^^ increases with the mass influx rate. 

Now, Bll assumes that M oc hB^^^; consequently, an argument similar to that above 
leads him to conclude that the mass pile-up at any such location is a runaway process: a 
net influx of material reduces hB^^^, which in turn yields (with his formula) a smaller M 
out of this location and hence even higher net influx, and so on. The general formula we 
use for the accretion rate, however, paints a more complex picture. Since /i-B^ax (and hence 
also r'^hB^^^, for a given r) decreases faster for higher influx rates, we have in general: 
d{r^hB^g^^)/dt = —f{dM/dr), where / is some positive function that captures the details 
of this evolution. With M oc d{r'^hB'^^^ / dr, we finally get: dMjdt cx —d[f{dM/dr)\/dr. 
If we ignore the outer differential, as Bll essentially does, then indeed M diminishes with 
time {dM/dt < 0) for a net mass influx {dM/dr > 0), leading to a runaway pile-up. The 
general equation, though, shows that whether M is enhanced or reduced with time depends 
on the radial variation of the mass influx rate dM/dr; the process only runs away if the latter 
increases with r, which is not necessarily true even if dM/dr > 0. Analogous arguments can 
be made for the slightly different phenomenologies that apply when the active layer extends 
to the midplane, and/or there is a net mass outflux instead of influx. In summary, the 
final outcome may be runaway mass pile-up in some regions, runaway evacuation in others, 
and oscillations about (or towards) a steady-state in yet others, depending on the precise 
variations (both first and second-order) of the accretion rate with radius (and the specifics 
of the conditions in the active layer, which determine the function /). Further investigation 
of this requires detailed time-dependent simulations, beyond the scope of this paper. 

(iv) Whatever the final, or time-dependent, structure of the disk might be, the above dis- 
cussion shows that it is almost certain to comprise regions of both over- and underdensity; 



•^Note that a quasi-steady ionization and chemical equihbrium is always expected, because the recombi- 
nation and field growth/diffusion timescales are much shorter than the secular timescales over which mass 
accretion occurs; see Bll. 
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achieving steady-state with a surface density that decreases monotonically with radius ap- 
pears highly unhkely. On the one hand, this imphes that the formation of such disk struc- 
tures, currently attributed to either sculpting by planets or photoevaporation, may be a 
natural outcome of the MRI-driven accretion itself. On the other hand, this may enhance 
the former two processes: photoevaporation will be more efficient in rarefied (more optically 
thin) regions, while grain growth can be accelerated in overdense regions (due to both the 
increased density, as well as heightened trapping of grains due to the pressure gradient re- 
versal caused by the local density enhancement, similar to the situation envisaged by Ida 
& Lin (2008) near the ice-line). Finally, both disk gaps and pressure gradient reversals due 



to overdensities can halt planetary migration (Alexander & Pascucci 2012 Ida & Lin 2008) 



possibly making Type I migration timescales long enough to be consistent with observations. 

(v) In general, our inferred accretion rates increase outwards with radius, before falling off 
again in the outer disk, broadly consistent with the results of Bll (except that our M often 
reverses direction in the outer disk as well). The mean empirical rate of infall onto young 
stars - ~ lO~^M0yr~^ for solar type stars and 10~^°MQyr~^ for very low mass ones (e.g., 
Mohanty et al. 2005 MuzeroUe et al. 2005 ) - is achieved in the innermost parts of our disks 
for dust-to-gas mass ratios of < 10~^ for 0.7 and < 10~^ for 0.1 respectively. Thus 
a factor of ~10-1000 depletion in the mass sequestered in well-mixed small grains, relative 
to the standard ISM ratio of 10~^, appears necessary to reproduce the observed stellar M. 
Given that grain growth up to millimeter sizes seems common in these stars even by ~1 Myr 
( [Andrews fc Wilhams|r2005| |2007t [Ricci et al.||2010a||b| [Mohanty et~aI]|2012D , this appears 
plausible: since the mass of a grain rises as the cube of its radius, very few large grains 
are necessary to significantly depress the mass locked up in small ones (e.g., for a 1000- fold 
depletion in the total O.l/im grain mass, only one in a billion grains in the depleted disk 
needs to be millimeter-sized). These levels of depletion in small grains are also consistent 
with mid-infrared observations of T Tauri disks (e.g. 



Furlan et al. 2006) 



In this context, we point out that the M we derive using Bll's formula are 10-20 times 
larger than he finds, for the case of a solar-type star (0.7 Mq) with 0.1 /xm grains and a 
dust-to-gas ratio of 10~^; for the grain-free case, our values with his formula are up to a 
100 times larger than his. In both cases, this is because our -Bmax are larger, because of 
higher electron abundances. In the first instance, this is due to our more detailed treatment 
of the X-ray fiux and Compton scattering, as well as some differences in the chemistry and 
surface density (see discussion in §6.2.1), leading to 3-5x larger Xe (and thus -Bmax)- In the 
grain-free discussed at length in §6.1, our Xg are considerably higher, due to the 

difference in our chemical networks. 



(vi) Finally, in spite of our somewhat larger -Bmax in the presence of grains, compared to the 
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results of Bll, our field strengths at 50-100 AU are still <10 mG, for 0.7 Mq with dust-to-gas 
ratios 10~^-10~^ (and lower yet for 0.1 Mq). We therefore agree with Bll that the current 
non-detection of polarized emission from grains aligned with the field is consistent with 
active MRI. Specifically, SMA observations by Hughes et al. (2009), with angular resolution 
corresponding to 50-100 AU, have failed to detect such emission at mm-wavelengths from 
two solar-type classical T Tauri stars. Their models suggest that most of this emission is 
expected from 10-100 /im grains, for which the minimum field strength for grain alignment 
is 10-100 mG. Our results, however, show that fields this strong would disrupt the MRI in 
the outer disk in the presence of even a tiny fraction of very small grains; the non-detection, 
implying weaker fields, is thus compatible with active MRI. 



7. Conclusions and Future Work 

We have presented numerical simulations of X-ray driven MRI in disks around both 
solar-type and very low mass stars, including the the effects of ambipolar diffusion, in the 
presence of 0.1 yum grains. In his seminal study, Bll undertook an analogous investigation 
of solar-type stars. In addition to extending his work to the very low mass stellar domain, 
we have incorporated several improvements to his treatment: a more sophisticated analysis 
of the X-ray ionization; a generalized formulation of the accretion rate; an investigation 
of various degrees of grain depletion between the limiting cases of a standard ISM dust 
abundance and no dust at all; and an examination of the active layer properties over the 
entire disk from 0.1 to 100 AU. Our choice of the field strength at every radius is governed by 
the assumption that the MRI generates a magnetic field of the strength that maximizes the 
height-integrated accretion stress, which is very close to the maximum possible field strength 
that still permits active MRI (the same choice as Bll makes, but for a different physical 
reason). Our main conclusions can be summarized thus: 

• Ambipolar diffusion is the primary factor determining the thickness and location of the 
MRI-active layer around both solar-type and very low mass stars, and quenches the MRI 
over most of the disk: the active region constitutes only 5-10% of the total disk mass in 
all cases. Previous suggestions that reduced surface density in low-mass disks (around less 
massive stars) may cause the entire disk to become active (e.g., Hartmann et al. 2006[ ) are 



invalid, precisely because ambipolar diffusion becomes stronger at lower densities. 

• Additionally, small well-mixed grains strongly suppress the MRI. Achieving the empirical 
accretion rates in solar-type and very low mass stars requires that the total mass in such 
grains be depleted (via grain growth and/or settling) by a factor of 10-1000, relative to the 
standard dust-to-gas mass ratio of 10~^. This degree of depletion is not implausible, given 
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the rapid grain growth observed in classical T Tauri stars, and is consistent with the range 
of depletion suggested by mid- infrared observations of these disks (Furlan et al. 2006[ ). 



• The accretion rates are radially variable, and may even change sign (direction) at some 
locations. This can produce both runaway pile-up or evacuation of material at some radii, 
and oscillations about a steady-state in others. Radially, the disk is thus expected to comprise 
regions of both over- and underdensity. These can accelerate photoevaporation as well as 
grain coagulation and planetesimal formation, and reduce the efficiency of Type I migration. 

• The current non-detection of polarized emission from field-aligned grains in the outer disk 
is consistent with active MRI in these regions. 

Various questions still remain open, however. Some of the most critical, offering prof- 
itable avenues for future investigation, are as follows. 

• Chemistry: As discussed above, we find (in agreement with Bll and PCll) that a very 
substantial depletion of small grains is required for X-ray driven MRI to reproduce the 
observed accretion rates in classical T Tauri stars. However, it is precisely at such low grain 
abundances that the derived MRI accretion rate becomes sensitive to the particular chemical 
network adopted (as indicated by the significant disparity between our M and that of Bll, 
in the dust-free case). Such severe small-grain depletion may also favour the dominance 
of atomic hydrogen over molecular, contrary to our and Ell's assumption (as discussed in 
§4.2). Thus, a more detailed investigation of chemical networks and hydrogen recombination 
channels at very low grain abundances is essential to confirm (or refute) the importance of 
X-ray driven MRI for disk accretion. 

• PAHs: In an elegant semi-analytic study, PCll showed that PAHs - macromolecular sized 
dust particles - may pose an even graver threat to X-ray-driven MRI than the small grains 
considered here. This is because the efficiency with which dust grains soak up electrons 
(and hence also the rate of ion recombinations on grain surfaces) increases with the total 
grain surface area, and for a fixed dust mass, the total area rises in inverse proportion to the 
grain radius. PCll consider PAHs of radius 6 A, which would enhance the area available for 
adsorption and recombination by a factor of 200 compared to our 0.1 /im grains, and thus, 
prima facie, enormously suppress the MRI. However, PCll neglect the competing fact that 
grains are also better tied to field lines as they become smaller: not because of a reduced 
grain mass (which is irrelevant as long as it far exceeds that of a neutral particle), but because 
of a reduced collision frequency with neutrals. Specifically, from equation (19), the ratio of 
the Hall parameters for PAHs and ions (a measure of their relative attachments to the field) 
is /Spah/A ~ (I-^pahI < CT'W >i)/i\Zi\ < (TV >pah)- At a fiducial disk temperature of ~300K, 
corresponding to around 1 AU for a solar-type star, equation (26) implies that the coUisional 
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rate coefficient for 6 A PAH grains becomes < av >pah ~ 10~^ cm^ s~^, comparable to that 
for ions (equation (25), for both molecular and metallic ions). Thus, even if most small 
PAHs carry a single charge (as PCll find), they are still just as bound to the field as 
the ions (which are predominantly singly- charged as well), at this temperature. In cooler 
parts of the disk, small PAHs are even better tied to the field than ions are (since the rate 
coefficient for grains goes as a/T while that for ions is constant); the same is true if the 
PAHs are multiply-charged. In this situation, PAHs complement ions in tying neutrals to 
the field through collisions, thereby reducing ambipolar diffusion. If this effect overcomes the 
increased recombination these particles facilitate, then PAHs can promote the MRI instead 
of hindering it. This very promising volte-face by PAHs richly deserves further study. 



FUV-driven MRI: Following the work of PCll, Perez-Becker & Chiang (2011b) showed 



that the MRI can thrive even in the presence of PAHs if the ionization is driven by stellar 
FUV photons instead of X-rays: while the FUV radiation does not penetrate as far into the 
disk, it ionizes carbon and sulphur atoms to yield electron and ion fractional abundances 
orders of magnitude higher than achieved with X-rays, completely dwarfing PAH effects. At 
radii > 10 AU, the maximum inferred M in this case is comparable to the mean observed 
rate of 10~^ Mq yr~^ for solar-type stars, though ambipolar diffusion again causes the rate 
to fall well below the observed mean in the inner disk regions. Thus FUV ionization may 
explain the empirical accretion rates of classical T Tauri stars, at least beyond 10 AU, without 
having to invoke significant depletion of small grains. However, it is unclear how much of the 
FUV radiation actually reaches the disk surface: in accreting stars, these photons are mainly 
generated in the accretion shocks on the stellar surface, and should be strongly attenuated on 
their way to the disk by the surrounding infalling accretion columns. An analysis accounting 
for this effect is required to better evaluate the viability of FUV-driven MRI. 

• M variability and time- dependence: We find the accretion rates to be radially variable; 
this is a general feature of layered accretion (e.g., Gammie| ( 1996[ )). In addition, we see that 
this not only leads to a buildup of mass in the inner disk (e.g., Gammie 1996, Bll), but 



can cause either a pileup or evacuation of material at several locations, depending on how 
M changes with radius. While the importance of mass pileup for non-steady accretion and 
disk instabilities is recognized, the effects of accretion-related radial oscillations in density 



on planet formation and migration have been little studied (but see Ida & Lin (2008) for 
a discussion of viscosity-related overdensity near the ice-line, and [Alexander fc Pascucci 



(2012) for the effect of photoevaporation-driven gaps on migration). Doing so first requires 
describing these density variations in more detail, through time-dependent MRI accretion 
simulations in the future. Such dynamic simulations are also required to better treat the 
time-dependent effects of grain growth and settling (e.g.. Turner et al. 2010). 
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• non-ideal MHD: Finally our investigation, like that of Bll, is based on the latest numerical 
simulations of non-ideal MHD, which still have a number of important uncertainties. Two 
in particular are worth mentioning. The first is that our ambipolar diffusion criterion for the 
MRI, equation (4), originates in the unstratified calculations by Bai & Stone (2011); whether 
this continues to hold in the stratified domain remains to be seen. Second, we have ignored 
the effects of Hall diffusion on the MRI, guided by the non-linear simulations of Sano fc Stone 



(2002). As Wardle & Salmeron (2012) point out, however, these calculations do not probe 
the "deep" Hall regime, where \rjH\ > f]p > v\/VL. Wardle & Salmeron's /mear analysis, with 
vertical magnetic fields, suggests that Hall diffusion in this regime can change the column 
density of the active layer by an order of magnitude or more (either increasing or decreasing 
the column depending on whether the field is parallel or antiparallel to the rotation axis). 
Non-linear simulations are now essential to unravel the true role of Hall diffusivity here. The 
bottom line is that analyses such as ours depend perforce on the results of non-ideal MHD 
simulations, and advances in the latter are vital to further improve our understanding of the 
viability and importance of the MRI for protoplanetary disk accretion. 
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Fig. 1. — Ionization rate {black dashed line) and fractional abundances of the various species 
in our chemical network {colored lines, specified in key in top left panel), as a function of 
column density at 1 AU, for various dust-to-gas ratios and two stellar masses. Top: 0.7 M©. 
Bottom: 0.1 Mq. Left to Right: dust-to-gas ratio = lO^^ 10"^ 10-^ 0. See §6.1. 



Fig. 2.— Same as Fig. 1, but at 10 AU. 
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Fig. 3. — Ionization rate ( {black line) and fractional electron abundance Xe {red: with dust- 
to-gas ratio = 10~^, blue: no-dust), at 1 AU for 0.7 M©, as a function of number density {top 
panel) and height {bottom panel). The dashed vertical line shows ( = 10""^'^. See §6.1. 
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Fig. 4. — Various field-dependent quantities, for 0.7 Mq. Top plot (a): At 1 AU. Bottom plot 
(b): At 10 AU. Within eacli plot, we show:- Top left: hB^ as a function of field strength 
B. The dashed vertical line marks the field strength at which hB^ is maximum; note that 
this is very close to the maximum B allowed (which is the end of the curve). We set Bmax 
equal to the B at which hB"^ is maximized. Top right: MRI-active layer {grey region) as 
a function of field strength and vertical height. A^ > 1 everywhere to the right of the red 
curve] h.K > 10 everywhere to the right of the grey vertical line] and ji/ ^uAn > 1 everywhere 
to the left of the blue curve and right of the aqua curve. The horizontal dashed line marks 
B = -Bmax! note that at this field strength, both the bottom and top of the active region are 
set by the ambipolar condition (blue and aqua curves). Bottom left: A_4 (red), Ak (grey) and 
/3//3min (blue) as a function of height, for B = i?max- A_4 > 1 and /3//3mm > 1 above the black 
dashed line, and A/c > 10 above the grey dashed line. Bottom right: Various resistivities as 
a function of height, for B = -Bmax- Vo in "i^^d, \riH\ in green {solid and dashed: positive and 
negative rjn respectively), tja in blue. See §6.2.1. 
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Fig. 5. — Various resistivity and activity zones over the entire disk, for 0.7 Mq and dust-to- 
gas ratio = 10^^. Top panel: Dominant resistivities (colors specified in plot key). Middle 
panel: Same, with Dead Zone {brown) and Undead Zone {grey) overplotted. Bottom panel: 
Same, with Zombie Zone {yellow) overplotted as well. Where the Zombie Zone overlaps with 
the dead and Undead Zones, the latter are shown as a picket-fence. The MRI-active zone 
lies within the black-ringed region. Note that the bottom of the active zone is always set by 
the Zombie Zone (ambipolar condition), in this and all following plots (not obvious at the 
smallest radii at this plotting resolution, but nevertheless true). See §6.2.2. 
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Fig. 6. — Same as Fig. 5, but for 0.7 and dust-to-gas ratio = 10 
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Fig. 8. — Same as Fig. 5, but for 0.7 and dust-to- gas ratio = (no dust). 
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Fig. 9. — Same as Fig. 5, but for 0.1 and dust-to-gas ratio = 10 
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Fig. 10. — Same as Fig. 5, but for 0.1 Mq and dust-to-gas ratio = 10 
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Fig. 11. — Same as Fig. 5, but for 0.1 and dust-to-gas ratio = 10 
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Fig. 12.— 



Same as Fig. 5, but for 0.1 Mq and dust-to-gas ratio = (no dust). 
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Fig. 13. — Various quantities relevant to the accretion rate, plotted as a function of radius, for 
dust-to-gas ratio = 10~^. Top panels: 0.7 Mq. Bottom panels: 0.1 M©. Left to right: -Bmax! 
^-^max) which sets the accretion rate (modulo fl) using the formula of Bai (2011); r^hB"^^^, 
the slope of which sets the accretion rate (modulo rQ) using our generalized formula; the 
accretion rate M. The black curve shows M with Bai's formula (always +ve), and the colored 
curves show M with our formula {red: positive (inward) M, blue: negative (outward) M). 
Our M are overwhelmingly -|-ve for this dust-to-gas ratio, but not for others (following plots). 
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Fig. 16. — Same as Fig. 13, but with dust-to-gas ratio = (no dust). 



